%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                          Transitional Dynamics                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ---------------------------- Description -------------------------------- 

    % This scripts uses the equilibrium conditions of the model to compute
    % the transition of the economy from the calibrated BGP to the
    % counterfactual BGP.

%--------------------------------------------------------------------------
%                            Preliminaries                                %
%--------------------------------------------------------------------------

    % Determine the number of periods after the shock takes place
    TimePeriods = p.Tperiods;   

%--------------------------------------------------------------------------
%                                 BGP                                     %
%--------------------------------------------------------------------------

    % Before we compute the transitional dynamics implied by the population
    % growth path, we need to compute some objects on the BGP that will be
    % used for calculation.

    % ---------------------- Variable paths -------------------------------

    load InfoMatrix.mat

    % Path for the inverse of the variety intensity - starts on the BGP
    LFN_t = zeros(1, TimePeriods);
    LFN_t(1) = LFN;

    % Path for the variety growth rate path (g^N_t) - starts on the BGP
    gN_vector = zeros(1, TimePeriods);
    gN_vector(1) = p.eta;

    % Initial guess for the production share path (s_t^P)
    prod_share_vector = InfoMatrix(4, :);

    % Path for the misallocation term (M_t) - starts on the BGP
    M_vector = linspace(M, M*0.995, TimePeriods);
    M_vector = max(M_vector, M_counter);

    % Path for the labor wedge term (Lambda_t) - starts on the BGP
    Lam_vector = linspace(Lam, Lam*0.95, TimePeriods);
    Lam_vector = max(Lam_vector, Lam_counter);

    % Path for the labor wedge term growth rate (g_Lambda) - starts on the BGP
    gLam = zeros(1, TimePeriods);
    gLam(1) = 0;

    % --------------------- Function b_t(Delta) ---------------------------

    % For an overview of the function b_t(Delta) on the BGP, see Section (RP-3.2.1)

    % ---------------------------------------------------------------------

    % Determine number of points in the productivity grid
    u_Delta_pts = 100;

    % Compute step size to have CES mark-up and lambda on the grid
    u_Delta_step = (p.s/(p.s-1) - p.l)/(u_Delta_pts-1);

    % Construct the productivity grid
    u_Deltagrid = -u_Delta_step * ones(1,u_Delta_pts);
    u_Deltagrid(1) = p.s/(p.s-1);
    u_Deltagrid = fliplr(cumsum(u_Deltagrid));

    % Preallocate a matrix to store this function along the transition
    u_func = zeros(u_Delta_pts, TimePeriods);

    % Obtain the function b on the BGP - see Section (RP-3.2.1)
    ConstC_t = ConstC*p.yearlength;
    In_t = p.In*p.yearlength;
    ODE_const_t = (1/(ConstC_t*(p.s-1))-p.s/(ConstC_t*(p.s-1)+In_t*(p.s-1)^2)+1/(ConstC_t+In_t*p.s))*(p.s/(p.s-1))^(-p.s-ConstC_t/In_t);
    u0_func = @(d) (min(d,p.s/(p.s-1)).^(-p.s).*(min(d,p.s/(p.s-1))/(ConstC_t+In_t*(p.s-1))-1/(ConstC_t+In_t*p.s)) + ODE_const_t*min(d,p.s/(p.s-1)).^(ConstC_t/In_t))/(M^(p.s-1)*Lam^(p.s));

    % Allocate the BGP function to the initial column of matrix
    u_func(:, 1) = u0_func(u_Deltagrid);

%--------------------------------------------------------------------------
%                           Transitional Dynamics                         %
%--------------------------------------------------------------------------

    % Using the initial guess provided above, we update the current
    % production share as a function.  We continue as long as the
    % difference between two iterations is large.

    % ---------------------------------------------------------------------

    fprintf('------------------------- Transitional Dynamics started ---------------------\n')

    % Create a matrix with all the information the function needs
    InfoMatrix = [eta_vector; LFN_t; gN_vector; prod_share_vector; M_vector; Lam_vector; gLam];

    % If the switch is turned on...
    if Compute_TD == 1

        % ... the code runs and saves Transitional Dynamics
        [InfoMatrix, u_func] = Trans_FP(p, InfoMatrix, u_func, u_Deltagrid, F1, FNC, Deltagrid);
        %save InfoMatrix InfoMatrix u_func

    else

        % ... the code loads the results already saved
        load InfoMatrix.mat

    end

    fprintf('-------------------------- Transitional Dynamics ended ----------------------\n')
    fprintf('\n\n')
